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Report  Number  8:  THE  THEORY  OF  STANDARDIZED  TIME  SERIES; 

by  Peter  W.  Glynn  and  Donald  L.  Iglehart 

Ch  e  ;  'V  -  . r  r  '•  - 

^  Building  on  work  of  SGHRUBEN,  we 1  have  developed  a  general  framework  for 
the  analysis  of  standardized  time  series.  Under  mild  assumptions  on  the 
output  process  (see  (3.1)),  the  method  of  standardized  time  series  produces 
asymptotically  valid  confidence  lntevals  for  steadv-state  parameters. 
However,  these  Intervals  are  asymptotically  larger  (see  (5.16))  and  more 
variable  (see  (5.33))  then  those  steady-state  Intervals  obtained  by  a  method 
which  consistently  estimates  the  appropriate  steady-state  variance  constant 
(such  as  the  regenerative  method).  In  this  sense,  standardized  time  series 
confidence  Intervals  are  asymptotically  less  desirable  then  those 
constructed  by  consistent  estimation. 

These  results  do  not,  however,  preclude  the  possibility  that  standard¬ 
ized  times  series  may  be  superior  in  certain  small  sample  context;  this 
remains  an  area  for  future  work.  7  ,  ■  >  , 
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INTRODUCTION 


Let  T  ■  (Y(t)  :  t  ^  0}  be  a  real-valued  stochastic  process 
representing  the  output  of  a  simulation.  To  incorporate  stochastic 
sequences  {Y^  :  n  ^  0}  into  our  framework,  we  set  Y(t)  ■  Yjtj»  where 
[t]  is  the  greatest  Integer  less  than  or  equal  to  t.  Frequently,  a 
simulator  is  interested  in  estimating  steady-state  parameters  associated 
with  T.  Recently,  SCHRUBEN  (1983)  proposed  a  new  class  of  procedures, 
based  on  standardized  time  series,  for  dealing  with  the  steady-state 
simulation  problem.  Our  goal,  in  this  paper  is  to  generalize  the  method  of 
standardized  time  series  and  to  study  the  structure  of  such  procedures. 
Section  2  reviews  the  basic  concepts  of  weak  convergence  upon  which  the 
method  of  standardized  time  series  is  based. 

In  Section  3  the  method  of  standardized  time  series  is  introduced  and 
its  basic  properties  are  Investigated.  Section  4  gives  examples  of 
standardized  times  series,  while  Section  5  discusses  the  asymptotic 
behavior  of  the  method.  Section  6  provides  a  short  summary  of  the  major 
results  of  this  paper. 

2.  HEAR  CONVERGENCE  OF  STOCHASTIC  PROCESSES 

Let  X  ■  {Xn  :  n  ^  1}  be  a  sequence  of  real-valued  random  variables 
(RV's).  The  sequence  X  is  said  to  converge  weakly  to  a  r.v.  X 
(written  X  X  as  n  *  «)  if 


(2.1) 


P{X  <  x)  P{X  <  x}  , 


as  n  ®,  for  every  x  which  is  a  continuity  point  of  P{X  •}.  This 
mode  of  convergence  is  frequently  used  in  the  study  of  simulation  output 
analysis  algorithms.  For  example,  the  central  limit  theorem  (CLT)  is  a 
weak  convergence  statement  concerning  a  sequence  of  normalized  partial 
sums. 

It  is  well  known  (see,  for  example.  Theorem  4.4.2  of  CHUNG  (1974)) 
that  the  requirement  (2.1)  is  equivalent  to  demanding  that 

(2.2)  Ef (X  )  -  Ef ( X)  , 

n 

as  n  +  <*>,  for  every  bounded,  continuous  function  f  :  R  *  R. 

The  method  of  standardized  time  series  requires  that  one  study  weak 
convergence  properties  of  random  elements  X^  corresponding  to  stochastic 
processes.  Since  a  stochastic  process  may  be  regarded  as  a  random 
function,  it  is  natural  to  generalize  to  the  case  where  the  X^'s  take 
values  in  a  function  space.  The  precise  space  that  we  shall  require  is 
C [ 0 , 1 J ,  the  elements  of  which  are  continuous  functions  x  :  [0,1]  +  R; 
see  p.  54-61  of  BILLINGSLEY  (1968)  for  a  thorough  description  of  this 
space. 

As  (2.2)  indicates,  the  notion  of  weak  convergence  depends  upon 
defining  a  suitable  class  of  continuous  functions.  For  x,y  e  C[0,1],  let 

p(x,y)  *  sup{ | x( t)  -  y ( t ) |  :  0  <  t  £  1}  ; 
p(x,y)  measures  the  distance  between  the  two  elements  x,y. 
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(2.3)  DEFINITION.  A  function  f  :  C(0,lj  -*■  R  is  said  to  be  continuous 


if  f(x  )  f(x)  as  n  -►  ®,  whenever  p(x  ,  x)  0  as  n  -*•  ®,  where  the 
n  n 

x  's  and  x  are  elements  of  C[0,1]. 
n 


By  analogy  with  (2.2),  we  can  now  define  a  notion  of  weak  convergence 

on  C[0,1].  Let  {X^  :  n  ^  1}  be  a  sequence  of  random  elements  taking 

values  in  C[0,1]  (in  other  words,  the  X  's  correspond  to  stochastic 

n 

processes  with  sample  paths  in  C[0,1J).  If  X  is  a  random  element  of 

C [ 0 , 1 J ,  then  {X  :  n  >  1}  is  said  to  converge  weakly  to  X  (written 
n  — 

X  X)  if 
n 

(2.4)  Ef (X  )  -*>  Ef ( X) 

n 


as  n  ♦  ®,  for  every  bounded,  continuous  function  f  :  C [ 0 , 1 ]  >  R. 

The  method  of  standardized  time  series  is  based  on  the  following 
important  result,  known  as  the  continuous  mapping  theorem  (CMT).  For  a 
(measurable)  function  h  :  C[0,1]  >  R,  let  D(h)  be  the  set  of  elements 
x  e  C ( 0 , I ]  at  which  h  is  discontinuous  (in  other  words,  x  <  D(h)  if 


there  exists  a  sequence  {x  }  e  C{0,1)  for  which  p(x  ,  x)  -*•  0  as  n  ♦  ®, 

n  n 


with  f(x  )  i  f(x)). 
n 


(2.5  )  PROPOSITION.  Suppose  X  ,  X  are  random  elements  of  C [ 0 , 1 ]  such 

n 

that  X  ^  X  as  n  -*•  ®.  If  P{X  e  D(h)}  ■  0,  then  h(X  )  h(X)  as 
n  n 


See  p.  31  of  [1]  for  a  proof  of  this  result.  Loosely  speaking, 

Proposition  2.5  says  that  if  can  be  weakly  approximated  by  X,  then 

the  real-valued  random  variable  h(X  )  can  be  weakly  approximated,  in 

n 

distribution,  by  h(X),  provided  that  h  is  suitably  continuous.  To  be  of 
practical  benefit,  of  course,  it  is  also  necessary  to  choose  h's  for 
which  the  distribution  of  h(X)  is  known. 

3.  STANDARDIZED  TIME  SERIES 

Let  Y  *  (Y(t)  :  t  >  0}  be  a  real-valued  (m -asurable)  stochastic 
process  representing  the  output  of  a  simulation.  To  apply  the  method  of 
standardized  time  series  to  the  output  process  Y,  it  is  necessary  to  make 
the  following  assumption: 

(3.1)  There  exist  finite  constants  p  and  a  (a  positive)  such  that 

X  oB 
n 

as  n  •»  ®,  where  B  is  a  standard  Brownian  motion,  and 

X  (t)  -  nl/2(Y  (t)  -  ut) 
n  n 

with 

nt 

Y  (t)  -  /  -Y(s)ds/n  ,  for  0  <  t  <  1  . 
n  0  “  ~ 

Note  that  X  and  B  are  both  processes  whose  sample  paths  lie  in 
n 

C [ 0 , 1 J ,  so  that  the  weak  convergence  required  by  (3.1)  is  assumed  to  take 
place  in  the  function  space  C ( 0 , 1 J .  A  variety  of  different  output 
processes  satisfy  (3.1). 


(3.2)  PROPOSITION.  Let  Y  be  a  stationary  (measurable)  ^-mixing 
process  (see  p.  178  of  [lj  for  a  definition)  satisfying: 

(i)  EY2(0)  <  « 

00  1/2 

(ii)  /  *  <t)dt<- 

0 

CD 

(iii)  /  cov(Y(0) ,  Y(t) )dt  >  0  . 

0 

2 

Then,  (3.1)  holds  with  p  =  EY(0)  and  a  =*  2  f  cov(Y(0),  Y(t))dt. 

0 

For  a  proof,  see  p.  178-179,  as  well  as  a  remark  on  p.  150,  of  [1] 
(for  extensions  to  non-stationary  processes,  see  p.  179-182  of  [1]). 

(3.3)  PROPOSITION.  Let  Y(t)  -  Y.  , ,  where  (Y  :  n  >  0}  is  a 

[t]  n  - 

strongly  mixing,  strictly  stationary  sequence  (see  HALL  and  HEYDE  (1980) 
for  a  definition)  satisfying: 

(i)  e|y0|2+6<- 

(ii)  l  <z(n)6/(2+6)  <  • 
n*0 

CD 

(iii)  E(Y  -r)2  +  2  I  E(Yn-r)  (Y  -r)  >  0, 

u  k-1  u  K 

where  r  *  EY^;  5  >  0,  and  { ct( n)  :  n  0}  is  the  sequence  of  mixing 
constants.  Then,  (3.1)  holds  with  p  =  r  and  a2  =*  E(YqT)2  +  2 
E(Y0-r)(Yk-D. 


5 


For  a  proof,  see  p.  132  of  [8]. 


(3.4)  PROPOSITION.  Let  Y(t)  =*  Y.  .  ,  where  {Y  :  n  >  0}  is  an 

It  J  n  - 

associated  sequence  of  strictly  stationary  r.v. 's  (see  NEWMAN  and  WRIGHT 
(1981)  for  definitions)  satisfying: 

(i)  Yq  is  a  non-degenerate  r.v. 

CD 

(ii)  l  E(Y  -r)  (Y. -r)  <  ”,  where  r  *  EY.. 
k-1  0  k  0 

Then,  (3.1)  holds  with  p  «  r  and  a 2  =  E(Yg-r)2  +  2  ^  E( Y^-r) (Y^-r) . 

For  a  proof,  see  (10]. 

(3.5)  PROPOSITION.  Let  T  be  a  (possibly)  delayed  regenerative  process 

with  regeneration  times  0  T(0)  <  T(l)  <  T(2)  <  •••  .  Set 

T(n) 

Y  (f)  -  /  f ( Y(s) )ds 

"  T(n-l) 

x  -  T(n)  -  T(n-l) 
n 

and  assume  that: 

(i)  eC yl( |f |)2  +  X2)  <  - 

(ii)  E(YL(f)  -  rxL)2  >  0,  where  r  -  EY^O/Ex  . 

2  2 

Then,  (3.1)  holds  with  u  =  r  and  a  *  E(Y  (f)  -  rx  )  /Ex  . 


For  a  proof,  see  FREEDMAN  (1967)  (the  argument  given  there  for  the 
Markov  chain  case  easily  extends  to  the  general  regenerative  setting). 

Let  h  :  C [ 0 , 1 ]  >  R  be  the  mapping  defined  by  h(x)  *  x(l);  in  other 
words,  h  evaluates  x  at  the  point  t  *  1.  It  is  trivially  verified  that 
h(x^)  -*■  h(x)  whenever  p(xn>x)  +0  as  n  -*■  ®.  Hence,  Proposition  2.5 
implies  that 

Xn(l)  oB(  1)  , 
or 

(3.6)  n1/2(Y(n)  -  p)  cB(l) 

as  n  ♦  Application  of  a  standard  converging-together  argument  (see  p. 
93  of  CHUNG  [  4  ] )  yields 

(3.7)  PROPOSITION.  Assumption  (3.1)  guarantees  that 

Y(n)  =■■*  u 


as  n  >  <*> . 

Thus,  (3.1)  suffices  to  guarantee  that  the  steady-state  estimation 
problem  for  T  makes  sense;  p  is  the  steady-state  parameter  which  the 
simulator  wishes  to  estimate. 

Note  that  the  CLT  (3.6)  could  be  used  to  obtain  confidence  intervals 
for  p,  provided  that  a  were  known.  As  Schruben  points  out  in  [12],  the 
principle  underlying  standardized  time  series  is  to  "cancel  out"  the  a. 

The  cancellation  procedure  involves  choosing  a  function  g  from  the 
class  ;  JM  is  the  class  of  (measurable)  functions  g  :  C [ 0 , 1  ]  -*•  R 
such  that: 
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V 


,  ,1^.  I  >v  ,  n  1  V 


(3.8) 


(i)  g (ax)  =  ag(x)  for  a  >  0,  x  €  C[0,1], 

(ii)  g(x  -  pk)  =  g(x)  for  p  e  ]R  and  x  e  C[0,1],  where 


k(c)  *  t, 

(iii)  P{g(B)  >  0}  =  1, 

(iv)  P{B  €  D(g) }  =  0. 


(3.9)  THEOREM.  Suppose  that  g  e  31  .  Under  Assumption  (3.1), 


(3.10) 


Y  (1)  -  p 
n _ 

*<V 


=*> 


B(l) 

g(B) 


as  n  -*•  <*>  . 


PROOF.  Let  h  :  C(0,1]  ■*  R  be  the  mapping  defined  by  h(x)  =  x(l)/g(x) 
for  g(x)  f  0  (and  zero  elsewhere).  Assumptions  (3.8iii)  and  (3.8iv) 
allow  one  to  verify  that  P{oB  e  D(h)}  =  0.  Thus,  Proposition  2.5 
guarantees  that 


h( X  )  =»  h(orB) 
n 

as  n  -*■  <*>.  By  (3.8i),  h(oB)  =  B(l)/g(B)  (recall  that  a"  >  0). 
Furthermore , 


h(Xn) 


1/2  - 

n  (Y  (1)  -  4) 
_ n _ 

g(n1/,2(Yn  -  ku) 


V1}  -  * 

g(Yn  -  ku) 


Yn(D  -  P 


g(Yn) 
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For  the  second  step,  observe  that  the  symmetry  of  H  proves  that 
H(z(g;  1-6/2))  -  H(-z(g;  1-5/2))  =  1-6.  Set  b  =  z(g;  1-6/2).  Then,  for 


any  a  £  H,  (5.5)  yields 

H(<r+2b)  -  H(a)  _<  1-5  . 

Thus,  in  order  that  H(P)  -  H(  a)  =  1-6,  it  must  be  that  p-a  _>  2b;  proving 
our  assertion. 

We  turn  now  to  the  choice  of  g  £  .  Our  goal  is  to  find  g 

minimizing 

(5.6)  6(g)  *  Eg(B)  •  z(g;  1-6/2)  . 

Note  that  the  criterion  (5.6)  is  scale-invariant. 

(5.7)  LEMMA.  For  b  >  0,  <!)(bg)  *  6(g). 

PROOP.  Note  that  z(g;  1-5/2)  solves 

1  -  6/2  -  P { B(  1 )  <  z(g;  1-5/2)  •  g(B)} 

-  P { B(  1 )  <  -g  z(g;  1-5/2)  •  b  •  g(B)} 

-  P { B ( 1 )  <  z(bg;  1-5/2)  •  b  •  g(B) } 

so  that  the  continuity  and  strict  monotonicity  of  H  imply  that 

z(gb;  1-6/2)  *  -g  z(g;  1-5/2)  . 

Relation  (5.6)  then  yields  the  lemma. 


1 

■ 

1 


22 


Clearly,  it  is  desirable  to  obtain  confidence  intervals  with  as  small 
an  expected  length  as  possible.  From  Proposition  5.1,  it  seems  reasonable 
to  therefore  choose  a,  8.  and  g  e  J  such  that  Eg(B)  •  (0-a)  is 
minimized. 

(5.3)  PROPOSITION.  Suppose  g  e  JX.  Then,  for  a  100(1-6)%  confidence 
interval,  (3- a  is  minimized  by  choosing 

8  »  z(g;  1  -  6/2) 

a  =  -0 

where  z(g;  x)  solves  the  equation  H(z(g;  x))  =  P{B(l)/g(B)  <  z(g;  x) }  = 

(in  other  words,  the  confidence  interval  should  be  centered  at  Y  (1)). 

n 

PROOF.  We  proceed  in  two  steps.  First,  for  any  a  e  K  and  b,y  _>  0,  it 
is  easily  verified  that 

(5.4)  $((a+2b)y)  -  <$(ay)  <  $(by)  -  £<-  by)  . 

Integrating  both  sides  of  (5.4)  with  respect  to  G(dy)  and  using  (3.14), 
we  get 

(5.5)  H(a+2b)  -  H(a)  H(b)  -  H(-b)  . 

Furthermore,  the  symmetry  of  <?  and  (3.14)  implies  that  H  is  also 
symmetric,  in  the  sense  that  H(b)  -  H(0)  »  H(0)  -  H(-b)  for  b  0. 


(5.1)  PROPOSITION.  Assume  g  e  31,  and  that  (3.1)  holds. 

(a)  If  g  is  non-negative,  then 

1/2 

lim  n  EL  >  /Eg(B)  •  (B-ot)  . 

-  n  — 

n  *  <*> 

(b)  If  {g(X  )  :  n  .  1}  is  uniformly  integrable,  then 

n  — 

1/2 

lim  n  L  =  oEg(B)  •  (B_ct)  • 
n  a> 

Assumption  (3.1)  and  Proposition  2.5  guarantee  that  if  g  e  31,  then 

(5.2)  g(X  )  =»  ng(B)  , 

n 

as  n  ■*  If  g  is  non-negative,  then  Fatou's  lemma  can  be  applied  to 
(5.2)  to  conclude  that 

Eg(B)  <  Urn  Eg(Xn)  , 
n  ♦  ® 

proving  (a).  On  the  other  hand,  it  is  well-known  (see  CHUNG  (1974),  p.  96) 
that  uniform  integrability  implies  that 

Eg(B)  -  lim  Eg(Xn)  , 
n  ♦  • 


proving  (b). 


(4.10)  EXAMPLE.  Let  b  :  C[0,1]  -*•  R  be  defined  by 

b(x)  =  x(t*)/(t*(l-t*))1/2 

where  t*  *  inf{t  >  0  :  x(t*)  =  M*},  M*  =  max{x(t)  :  0  _<  t  <  1}.  SCHRUBEN 

2 

(1982)  showed  that  (b  °r)(B)  has  a  chi-square  distribution  with  3  degrees 
of  freedom.  Consequently, 


so  that 


-  -  3)  .  2  .1/2 

/m  B  '  (x3m) 


B  ( 1 )  3)  . 

-  -  t,  , 

✓m  g  (B)  3m 
ra 

where  t^  is  a  Student' s-t  RV  with  3m  degrees  of  freedom.  Confidence 

<*"N» 

intervals  based  on  g  (B)  as  defined  above  are  the  standardized  maximum 

m 

Intervals  of  [12]. 


5.  ASYMPTOTICS  FOR  STANDARDIZED  CONFIDENCE  INTERVALS 

In  this  section,  we  study  certain  asymptotic  properties  of  standard¬ 
ized  confidence  intervals.  In  particular,  we  consider  the  asymptotics  of 
the  expected  length  of  such  conf-idence  intervals,  as  well  as  the  end-point 
variability  of  these  intervals. 

Now,  from  (3.15),  it  is  clear  that  the  width  of  the  interval  (3.15)  is 


given  by 


/12m  g  (B)  -  (  l  bJ(1))  7  =  (y2)1/2 

m  i-0  ! 


where  x  denotes  a  chi-square  RV  with  m  degrees  of  freedom.  The 
m 

<v 

chi-square  property  of  g  (B)  makes  standardized  time  series  based  on 

m 

g  particularly  attractive,  since  in  that  case 
m 


B(  1) 


3) 


( 12m) 1/2  I  (B)  m 
m 


where  t  is  the  Student's-t  distribution  with  m  degrees  of  freedom;  the 
m 

limit  Theorem  (3.10)  can  then  be  used  to  construct  confidence  intervals  for 
H.  These  confidence  intervals,  which  were  suggested  by  SCHRUBEN  (1983), 

_  i/2 

are  based  on  the  so-called  standardized  sum  process  (Y  (1)  -  u)/(12m) 

n 

g(Y  (•))• 
m  n 


(4.9)  EXAMPLE.  The  map  b  :  C ( 0 , 1  ]  ♦  R  defined  by 


1 

b(x)  -  /  jx(t) j  dt 
0 


also  lies  in  the  class  37.  Furthermore,  the  distribution  of  (bO)(B) 

is  known;  see  JOHNSON  and  KILLEEN  (1983).  However,  the  distributions  of 

*  — 

both  g^B)  and  g^(B)  are  quite  complicated,  and  this  would  appear 
to  limit  the  applicability  of  this  method. 
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.1 


To  calculate  the  distribution  of  g^B)  and  gm(B),  convenient 

to  first  find  the  distribution  of  (b°D(B).  Note  that  the  continuity  of 
B  implies  that 

l  m  1 

(4.8)  ±  l  (TB) (i/m)  >  /  (TB)(t)dt 

m  i-1  0 

as  m  >  a.s.  The  left-hand  side  of  (4.8)  is  normally  distributed  with 
mean  zero  and  variance 


m  m 

-2  I  l  covf (rB)(k/m) ,  ( rB) ( A/m)  ] 
a  k=*l  A-l 


1  I  I  fmin(^,  ~)  *  v(m) 

m2  k-1  A-l  ®  °  m2  ’ 


Note  that  v(m)  is  a  Riemann  approximation  to  the  integral 


1  1 

/  /  [min(s.t)  -  st]ds  dt 

0  0 

which  has  value  1/12.  Thus,  v(m)  1/12;  hence,  taking  characteristic 
functions  of  both  sides  of  (4.8)  shows  that  the  right-hand  side  of  (4.8)  is 
normally  distributed  with  mean-zero  and  variance  1/12.  Since  A^B ,  .... 
A^jB  are  independent  Brownian  motions,  it  follows  that 


<2> 


m-1 


'12m  g  (B)  =-  [  |B  (1)  j 

i-0  1 


where  B  ,  ...,  B  are  independent  standard  Brownian  motions;  on  the 
u  m- 1 

other  hand. 
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intervals  converge  to  those  associated  with  a  normal  approximation.  This 
phenomenon  is  consistent  with  that  observed  in  the  method  of  batch  means, 
where  it  is  known  that  as  m  the  Student' s-t  distribution  approaches 

a  normal. 

A  second  "extension"  method  involves  defining  the  class 
«  (b  eUl  :  P(B  €  D(b2°T)}  -  0}  , 

2 

where  b  (x)  -  b(x)  •  b(x).  The  following  proposition  has  a  proof  similar 
to  that  of  Proposition  4.2. 

(4.5)  PROPOSITION.  If  b  £  37.,  then  g  e  ^W(m  >  1),  where 

2  m  — 

2  v  1/2 

g_  =  (  I  b  «r»A  ) 
m  i=0  1 

The  analogue  to  Proposition  4.4  Is  then  given  by 

(4.6)  PROPOSITION.  g  (B)  =*  (E( b2°r°B) ) 1/2  as  m  ♦ 

m 

Thus,  confidence  intervals  based  on  g  ( •)  will,  for  large  m, 

m 

correspond  to  that  associated  with  a  normal  approximation.  We  now  turn  to 

*  "*  , 

some  specific  examples  of  gm ' 3  and  s. 

(4.7)  EXAMPLE.  Let  b  :  C ( 0 , 1 J  *  R  be  defined  by 

1 

b(x)  *  |/  x(t)dt  J 
0 
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*  1/2 

(4.4)  PROPOSITION.  g  (B)/m  '  E(b°roB)  as  m  -*>  *. 

m 

PROOF.  Note  chat 

g*(B)/m1/2  -  -  l  b»r°(m1/2AiB)  . 
m  m  i-o  1 

1/2 

But  {m  A^B  :  0  i  <  m}  has  the  same  distribution  as  a  collection  of  m 
independent  standard  Brownian  motions  {B^  :  0  £  i  <  m},  so 

*,_w  1/2  3)  1  “v1  .  _  B 

g*(B)/m  =*  —  )  b°r°B 

m  i=0  1 

0®  denotes  equality  in  distribution;.  Of  course,  the  strong  law  of  large 
numbers  guarantees  that 

m-1 

-  I  b->r«B.  E( b»r°B)  a.s. 

■  1-0  1 

as  m  -*■  ®,  proving  the  result. 

Thus,  if  E( b°r<»B)  <  «,  we  observe  that 

■*/2».»>  -  ,  y'hm 

g*(Y  )/E<b°r°B)  g*(B)/E(b«roB) 

m  n  m 

as  n  -*•  ®,  where  the  limit  RV,  for  large  m,  is  the  normally  distributed 
quantity  B(l).  Hence,  as  the  method  of  standardized  time  series  is 
extended  to  more  and  more  increments,  the  corresponding  confidence 
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It  is  evident  from  (4.3)  that 


g  *  I  b°*  *  I  bo<?  or  , 

m  i-1  i-1 


so  that  if 


b  -  l  bo*  , 
■  i-I  1 


*  * 

we  have  a  representation  of  g  of  the  form  g  *  b  »F. 

m  mm 

Clearly,  b  satisfies  (3. lli).  For  (3.11ii),  observe  that  A  B  is  a 
m  i 

Brownian  motion  so  (3.11ii)  implies  that  (b*T®A^)(B)  >  0  a.s.  for 

0  <  i  <  m,  thus  yielding  (3.11ii)  for  b  . 

—  m 

For  (3.11iii),  note  that  the  continuity  of  A^  implies  that 


D(g  )C  U  (x  :  A  x  £  D(boT)} 
m  i-0  1 


so  that 


P{B  c  D(g  )}  <  [  P  {A ,  B  €  D(bor)  } 

m  ”  i-0  1 


But  A^B  is  a  Brownian  motion  so  that  (3. 1 liil)  shows  that 

P{A  B  c  D(boT)}  *»  0  for  0  <  i  <  m,  yielding  (3.  lliii)  for  g  . 
i  —  m 

* 

It  is  of  some  interest  to  consider  the  behavior  of  g  (B)  for  large  m 

m 
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As  we  have  already  seen,  the  fundamental  assumption  of  the  method  of 
standardized  time  series  Is  that  the  output  process  may  be  approximated  by 
a  Brownian  motion.  Intuitively,  then.  It  should  follows  that  the 
Increments  of  the  output  process  can  be  approximated  by  the  increments  of 
Brownian  motion.  This  suggests  that  one  might  try  to  extend  the  power  of 
the  method  of  standardized  time  series  by  applying  the  procedure  separately 
to  each  increment  of  the  output  process,  and  then  "patching"  the  Increments 
together.  In  some  sense,  this  phenomenon  occurs  in  the  method  of  batch 
means,  and  is  related  to  the  somewhat  arbitrary  nature  of  the  parameter 
m.  In  any  case,  we  now  present  two  “extension"  methods. 

Let  A^  :  C[0,1]  •*  C[0,1]  be  the  map  defined  by 

(Ajx)(t)  -  x((i+t)/m)  -  x(i/m)  ,  0  t  <  1 

for  0  <  1  <  m  (m  _>  1) ;  the  key  to  our  first  “extension"  procedure  Is  the 
following  result. 


(4.2)  PROPOSITION.  If  b  e  71,  then  g*  e  JW  (m  >  1),  where 

m  — 


m- 

-  I  b°r®A 
1=0 


PROOF.  We  shall  show  g  can  be  represented  as  g_  *  b  °r,  where 

m  Ttt  m 

b  e  Jl,  thereby  proving  that  g  e  JW*.  Let  V.  :  C [ 0 , 1 J  C [ 0 , 1  ]  be 
m  mi 

given  by 


(’?ix){t)  -  x((i+t)/m)  -  x(i/m)(l-t)  -  tx((i+l)/m) 


for  0  <  t  <  1  (0<^i<  m).  The  following  relations  are  easily  verified: 


g  (x)  *  f(-^r)  l  (A  x(i/m)  -  x(l)/m)2  11/2  , 

m  m- l  m 

where  A,  x(t)  ■  x(t)  -  x(t  -  1/h).  Note  that  A  B(i/m)  (i  *  1,  m) 

h  m 

are  increments  of  standard  Brownian  motion,  and  are  therefore  independent 

and  identically  distributed  normal  RV's  with  mean  zero  and  variance  1/m. 

Also,  B(l)/m  is  the  sample  mean  of  these  increments.  Hence  B(l)/g  (B) 

m 

has  a  Student' s-t  distribution  with  m-1  degrees  of  freedom. 

On  the  other  hand, 


where 


*.<V  ■  "'1/2t^r  I,  (Vn>  Z,  Vn>,2] 

i=l  j*l  J 


2 -,1/2 


in/m 

Z  (n)  =  /  Y(s)ds/(n/m) 

(i-l)n/m 


is  the  i'th  batch  mean  of  the  process  {Y(t)  :  0  £  t  <  n}.  Specializing 
Theorem  3.9  to  our  example  therefore  allows  us  to  conclude  that 


Zj  (n)  )2  ]1/2  »» 


as  n  *  ®,  where  t  .  is  a  Student' s-t  RV  with  m-1  degrees  of 

m—  1 

freedom.  To  summarize,  we  have  just  shown  that  the  method  of  batch  means, 
with  the  number  of  batches  fixed  at  m  2,  is  asymptotically  valid  under 
condition  (3.1).  This  result  complements  a  similar  theorem  due  to 
BRILLINGER  (1973). 
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P{(Yn(l)-n)/g(Yn)  <  x}  +  H(x) 


as  n  -*•  ®  for  all  x  e  R.  Hence,  to  obtain  a  100(1-6)%  confidence 
interval,  one  selects  a  and  (5  such  that  H(  p)  -  H(  a)  m  1-6  (such  or,  8 
exist  since  H(  •)  is  continuous;  also,  (3.14)  implies  that  H( •)  is 
strictly  increasing).  Then,  the  interval 

(3.15)  [Yn(l)  -  g(Yn)B,  Yn(l)  -  g(Yn)a] 

is  an  asymptotic  100(1-6)%  confidence  interval  for  p. 

The  process  (Y  -  pk)/g(Y  )  is  called  a  standardized  tine 
n  n 

series.  Theorem  3.9  and  Proposition  3.12  show  that  every  b  <z  71  gives 
rise  to  a  particular  standardized  time  series  procedure;  (3.15)  is  then  the 
corresponding  confidence  interval  for  p. 

4.  EXAMPLES  OF  STANDARDIZED  TIME  SERIES 

Our  first  example  of  a  standardized  time  series  captures  a  methodology 
which  has  been  extensively  studied  in  the  simulation  literature,  namely  the 
method  of  batch  means. 

(4.1)  EXAMPLE.  Let  b  :  C[0,1]  -*>  R  be  defined  by 

m 

b  (x)  -  r(-=r)  I  (x(i/m)  -  x((i-l)/m))2  11/2  , 

m  m-1  ±ml 

for  m  >  2.  It  is  easily  verified  that  b  e  Ul ,  so  that  g  *  b  °F  e 
—  m  m  m 

(see  Proposition  3.12).  But 


(T(x-pk))(t)  -  x(t)  -  pk(t)  -  t(x(l)  -  pk(l)) 


=*  x(t)  -  tx( i)  *  (rx)(t) 

so  T(x-pk)  ■  Tx;  hence,  g(x-pk)  ■  g(x). 

To  prove  that  SM  c_  JM* ,  consider  g  e  Jt .  We  claim  that  g  can  be 
represented  in  the  form  g  =  b°r,  by  setting  b  ■  g.  Recall  that  g(x)  * 
g(x-pk)  for  all  p  e  R.  In  particular,  setting  p  »  x(l),  we  see  that 
g(x)  *  g(Fx),  proving  our  assertion. 

We  can  now  obtain  the  following  result. 

(3.13)  PROPOSITION.  If  g  e  SH,  then  B(I)  Is  independent  of  g(B). 

PROOF.  It  is  well  known  that  the  process  B(t)  -  tB(l)  (0  _<  t  <  1)  is 

independent  of  B(l)  (see  p.  84  of  [1],  for  example).  In  other  words,  13 
is  independent  of  B(l),  which,  of  course,  Implies  that  g(B)  *  (b»r)(B) 
is  independent  of  B(l). 

Let  ■Kx)  -  P(B(  1)  <  x},  G(x)  -  P{g(B)  <  x},  and  H(x)  -  P{B(l)/g(B)  <  x}. 

Then, 

00 

(3.14)  H(x)  *  /  <5(xy)  G(dy) 

0 

by  Proposition  (3.13).  The  continuity  of  <5(  •)  and  the  bounded 
convergence  theorem  imply  that  the  right-hand  side  of  (3.14)  is  continuous 
everywhere  in  x.  Thus,  by  (2.1)  and  (3.10),  it  follows  that  under  the 
conditions  of  (3.1), 
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where  the  last  equality  is  due  to  (3.8ii).  These  observations  immediately 
yield  the  theorem. 

The  proof  clearly  indicates  the  role  of  Assumption  (3.8ii);  this 
condition  guarantees  that  g(X^)  does  not  depend  on  the  unknown  parameter 
M- 

To  construct  confidence  intervals  based  on  (3.10),  we  need  to  learn 
more  about  the  limit  RV  B(l)/g(B).  We  start  by  obtaining  an  alternative 
description  of  3t  .  Let  r  :  C 1 0 , 1 }  ■*•0(0,1]  be  the  map  define  by 

(Fx)(t)  -  x(t)  -  tx( 1 )  . 

Let  JI  be  the  class  of  functions  b  :  C[0,1]  +  F  which  satisfy: 

(3.11)  (i)  b(  coc)  *  ab(x)  for  a  >  0,  x  e  C[0,1], 

(ii)  P{(b<>r)(B)  >  0}  =  1, 

(iii)  P{B  e  D(b»r)}  -  0. 

Set  JM*  =  {g  :  g  ■  b»r,  b  e  72}. 

(3.12)  PROPOSITION.  31*  -  3t : 

PROOF.  We  first  show  that  JW*  =_  JW.  Suppose  that  g  -  b°r,  where 

b  e  71.  Clearly,  g  satisfies  (3.8i),  (3.8iii),  and  (3.8iv).  For  (3.8ii), 
observe  that 
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(5*8)  THEOREM.  Suppose  g  e  31 .  Then 


-1, 


*(g)  >  $  (1-6/2) 


where  <5  *  is  the  inverse  of  the  normal  cumulative  distribution  function 

«. 

PROOF.  By  Lemma  5.7,  we  may  scale  g  so  that 

z(g;  1-6/2)  »  1  . 


Now,  (5.9)  implies  that 


H( 1)  -  1-6/2  , 
or 

CO 

f  $(y)  G  (dy)  -  1-6/2  , 
0  g 


where  G  (dy)  ■  P{g(B)  c  dy}  (see  (3.12)).  Thus,  we  are  to  show  that 
8 


(5.10)  4>( g)  ■  Eg(B)  •  z(g;  1-6/2) 


/  (1-  G  (y))dy  >  $_1(l-6/2)  , 


subject  to 
(5.11) 


/  <Ky)  G  (dy)  -  1-6/2  . 
0  8 


Integrating  by  parts,  we  find  that 


23 


=■  /  G  (y)  <J>(y)dy  , 
0  8 


where 


Vy)  ’ 


Let  K(y) 


1  -  G  (y)  and  4>(y)  is  the  normal  density  function 

s 

be  the  distribution  function  defined  by 


0;  y  <  p 

K(y)  *  { 

1;  y  >.  P 

where  p  *  $  *( 1-6/2).  Note  that 


/  K(y)  4>(y)dy  -  1-6/2 
0 

and 


/  K(y)dy  -  p  , 
0 


where  K(y)  *  1  -  K(y).  Thus,  we  can  reformulate  (5.10)  and  (5.11) 
show  that 


(5.12)  /  (G  (y)  -  K(y) )dy  _>  0 

0  g. 

subject  to 

CD 

(5.13)  /  (G  (y)  -  K(y) )  <t>(y)dy  -  0  . 

0  8 


as 


Since  0  is  strictly  decreasing  on  [0,®),  and  because 


it  follows  that 


P  _  P  _ 

(5.14)  /  (G  (y)  -  K(y) )dy  •  6(p)  >  f  ( G  (y)  -  K(y))  <j)(y)dy  . 

0  8  0  8 


and 


(5.15)  /  (G  (y)  -  K(y))dy  •  4>(p)  >  f 

P  8  P 


Adding  (5.14)  and  (5.15)  together,  we  get 


/  (G  (y)  -  K(y))dy  •  <)>(p)  >  / 
0  8  0 

Relation  (5.13)  then  yields  (5.12). 


(G  (y)  -  K(y) )  4»(y)dy  . 

O 


(Gg(y)  -  K(y) )  <()(y)dy  . 


(5.16)  COROLLARY.  Suppose  g  e  31  is  non-negative.  Under  Assumption 
(3.1), 


1/2  -1 
11m  n  '  EL  >  2a  «  (1-6/2)  . 

n  +  « 

This  corollary  follows  immediately  from  Propositions  5.1  and  5.3,  and 
Theorem  5.8.  The  lower  bound  of  Corollary  5.16  has  an  important 
interpretation.  Consider  a  steady-state  simulation  output  analysis 
algorithm  which  is  based  on  constructing  an  estimator  s^  which 
consistently  estimates  o: 


as  n  -*■  ®.  Among  Che  algorichms  of  Chis  type  are  Che  regeneraCive  meChod 
of  simulaClon,  specCral  mechods,  and  auCoregressive  procedures  (see 
ChapCer  3  of  BRATLEY,  FOX  and  SCHRAGE  (1983)  for  a  descripCion  of  Chese 
Cechniques).  The  following  proposicion  is  a  scraighcforward  applicaCion  of 
Che  converging-CogeCher  lemma  (see  p.  25  of  [1]). 

(5.18)  PROPOSITION.  If  s  is  an  esCimaCor  saCisfying  (5.17),  chen 

n 

(3.1)  implies  ChaC 

(5.19)  n1/2(Y  (1)  -  p)/s  B ( 1 )  , 

n  n 

as  n  ■>  ®. 

The  weak  convergence  resulC  (5.19)  permiCs  consCrucCion  of  asympcocic 
100(1-6)X  confidence  incervals  for  p: 

s  s 

(5.20)  [Yn(l)  -  z(6)  -jjj,  Yn(l)  +  a(  6)  -fij]  , 

n  n 

where  z(6)  -  <5  *(1-6/2).  If  L  is  Che  lengCh  of  Che  iiCerval  (5.20),  ic 

n 

is  clear  chac  as  n  ♦  ®, 

(5.21)  n1/2  L  *+  2o  $“l(l-6/2)  , 

n 

which  is  precisely  Che  lower  bound  of  Corollary  5.16.  If  (s^;  n  1 }  is 
uniformly  inCegrable  (condiCions  guaranCeeing  Chis  appear  in  GLYNN  and 
IGLEHART  (1985b),  SecClon  6),  Chen  we  furCher  have  ChaC 


(5.22)  lim  n  '  EL  -  2a  $  (1-6/2)  . 

n  •>  »  ® 

Corollary  5.16,  and  the  limit  theorems  (5.21)  and  (5.22)  suggest  that, 
from  the  viewpoint  of  expected  confidence  interval  length,  output  analysis 
methods  which  consistently  estimate  a  dominate  standardized  time  series 
procedures  asymptotically. 

One  further  point,  pertinent  to  expected  confidence  interval  length, 

remains  to  be  investigated.  The  examples  of  Section  4  show  that  for  any 

k  1,  there  exists  g^  e  JM  such  that  B(l)/g^(B)  has  a  Student's-t 

distribution  with  k  degrees  of  freedom.  If  g  (X  )  is  uniformly 

k  n 

integrable,  then  it  follows  that  if  L^(k)  is  the  length  of  such  a 
confidence  interval, 


lim  n1^2  EL  (k)  *  2a  H,  1  ( 1—6/ 2)  , 
n  k 


where  H^(p)  is  the  pth  quantile  of  a  Student's-t  with  k  degrees  of 
freedom.  Since 


lim  H^(l-6/2)  -  ®"1(l-6/2)  , 
k  ® 


this  discussion  suggests  that 


1/2  -1 

(5.23)  inf  Urn  n  '  EL  -  2a  <5  (1-5/2)  ; 

g  £  J  n  +  ®  n 


27 


thus,  the  lower  bound  of  Corollary  5.16  is  tight.  Relation  (5.23)  raises 
the  question  of  whether  there  exists  g  e  SH  such  that 

1/2  -1 

(5.24)  lim  n  '  EL  -  la  $  (1-6/2)  ; 

-  n 

n  >  ® 

in  other  words,  is  the  lower  bound  attained  within  SH  ? 

A  glance  at  the  proof  of  Theorem  5.8  shows  that  1 

6(g)  >  $“1(l-6/2) 

unless  G  (dy)  is  a  point-mass  distribution.  Thus,  in  order  to  find 
8 

g  e  JJK  satisfying  (5.24),  it  oust  be  that 

(5.25)  P{g(oB)  *  aa}  ■  1 

for  some  a  >  0.  Our  next  result  shows  that  such  a  g  cannot  exist. 

(5.26)  PROPOSITION.  There  exists  no  g  e  31  such  that  (5.25)  holds. 

PROOF.  We  will  prove  something  stronger:  the  requirements 

P{ B  c  D(g) }  -  0  and  (5.25)  are  incompatible.  We  start  by  showing  that  for 

every  x  €  Cq [0,1]  =  {x  €  Cl0,l]  :  x(0)  ■  0}  and  e  >  0, 

(5.27)  P{p(oB,  x)  <  e}  >  0  . 

To  see  this,  fix  x  «  Cq[0,1]  and  6  >  0.  Since  [0,1]  is  compact,  x  is 
uniformly  continuous  on  [0,1],  so  there  exists  N  ■  N( e)  such  that 
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(5.28) 


|x(t)  -  x(k/N) |  <  e/4 


for  kN  £  t  <_  (k+l)/N,  where  0  £  k  <  N.  Now,  the  independent  increments 
of  Brownian  motion  imply  that  if  Az(k/N)  =  z((k+l)/N)  -  z(k/N) ,  then 


(5.29)  P(A( e))  =  P j | <jAB(k/N)  -  Ax(k/N)|  <  e/4N, 

a  max  |B(t)  -  B(k/N) j  <  e/2,  0  <  k  <  N} 

k/N<t<(k+l)/N 

N-l 

*  IT  P{  |  aAB(k/N)  -  Ax(k/N)|  <  e/4N, 
k-0 

a  max  { B( t )  -  B(k/N)|  <  e/2}  >  0  , 

k/N<t<(k+l)/N 


by  virtue  of  the  fact  that  for  any  z  with 

maX0<s<t  ^  >  °*  Now*  on  the  event 

inequality  argument  shows  that 


[z|  <  p,  P { J B( t)  -  z|  <  r i, 
A(  e) ,  a  simple  triangle 


}aB(t)  -  x(t) |  <  e 


for  0  _<  t  _<  1  (use  (5.28)),  proving  (5.27). 

From  (5.27)  and  (5.25),  it  follows  that  for  some  a  >  0 


P{p(aB,x)  <  e,  g(aB)  ■  aa}  >  0  , 

so  that  there  necessarily  exists  y  ■  y(x,e)  such  that  p(y,x)  <  e  with 
g(y)  *  an.  Thus,  the  range  of  g  over  any  e-neighborhood  of  x  contains 
the  set  {at a:  a  >  0};  clearly,  then  g  can  not  be  continuous  at  x. 


Hence,  x  €  D(g).  Since  x  was  arbitrary,  this  implies  that  D(g)  * 
Cq 10,1] ,  violating  the  assumption  P{B  e  D(g)}  ■  0. 

We  now  turn  to  the  question  of  end-point  variability.  To  be  precise, 
observe  that  if  g  €  JA,  then  (3.1)  implies  that 

(5.30)  n1/2  L  ^  ag(B)  •  ( 0-  a) 

as  n  «  (see  (5.2)  for  a  more  complete  argument).  The  limit  distribu¬ 
tion  of  the  confidence  interval  length  is,  of  course,  degenerate  if  (and 
only  if)  g(B)  is  degenerate.  Suppose  that,  in  fact,  g(B)  is  degenerate 
so  that  there  exists  a  such  that  P{g(B)  ■  a}  ■  1.  Note  that  a  >  0  by 
(3.8iii).  On  the  other  hand,  it  follows  from  (3.8i)  that 

(5.31)  P{g(oB)  »  aa}  »  1 

for  all  a  >  0.  But  (5.31)  is,  of  course,  just  (5.25);  Proposition  5.26 
therefore  proves  that  no  such  g  can  exist.  Consequently,  we  may  conclude 
that  g(B)  must  be  non-degenerate.  The  limit  theorem  (5.30)  therefore 
states  that  exhibits  non-degenerate  random  fluctuations  of  order 

n-1/2 

n- 

Another  way  to  quantify  the  above  phenomenon  is  to  examine  the 
? 

quantity  E(L  -  EL  )  . 

n  n 

(5.32)  PROPOSITION.  If  (g2(Xn):  n  >  1)  is  uniformly  integrable,  then 
under  (3.1), 


(5.33) 


lim  nE(L  -  EL  )2  -  c2E(g(B)  -  Eg(B))2  (g-a)2 
n  ®  n  n 

provided  g  e  31 .  Furthermore,  the  right-hand  side  of  (5.33)  is  positive. 

2 

PROOF.  The  uniform  integrability  of  {g  (X^):  n^  1}  implies  that  of 
(g(Xn):  n  1}  (see  p.  100  of  [4]),  so 

1/2 

lim  n  Eg(Yn)  *  lim  Eg(Xn)  -  Eg(B)  , 
n  +  ®  n  -*■  ® 

and 

lim  n  Eg2(Yn)  -  lim  Eg2(X  )  -  Eg2(B)  ; 
n  >  ®  n  “ 

combining  the  above  two  limit  relations  yields  (5.33).  As  for  the 
positivity,  this  follows  from  the  non-degeneracy  of  g(B)  for  g  e  31. 

We  now  wish  to  compare  the  end-point  variability  of  standardized  time 
series  procedures  to  that  obtained  via  methods  which  consistently  estimate 
o.  Our  analysis  will  be  restricted  to  the  regenerative  method  of 
simulation;  we  do  this  only  because  the  required  limit  theorems  are 
available  in  this  context. 

1/2 

As  (5.21)  indicates,  n  converges  to  a  degenerate  r.v.  Thus, 

-1/2 

L^  asymptotically  exhibits  no  random  fluctuations  of  order  n  .  We 
can,  in  fact,  be  more  precise. 

(5.34)  PROPOSITION.  Let  T  be  a  regenerative  process  satisfying 

E( (Yj | f | )^  +  t^)  <  ®  (see  (3.3)  for  the  definition  of  Y^( j  f | )  and 

t  ).  Then,  if  s  is  the  regenerative  estimator  for  a,  there  exists  n 
n  n 

such  that 
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as  n  -►  ®  ; 


(5.35)  (i)  n(L  -  EL  )  -»  nN(0,l) 

n  n 

(ii)  n2E(L  -  EL  )2  r?  as  n  -*•  ®  . 

n  n 

PROOF.  Under  the  above  moment  hypothesis,  there  exists  <  such  that 

1/2 

(5.36)  n  '  (s  -a)  =»  <N(0,1)  , 

n 

2 

as  n  *  <*>;  furthermore,  the  sequence  {n(sn~a)  :  n  ^  1}  is  uniformly 
integrable  (see  Sections  5  and  6  of  (7j).  Thus, 

(5.37)  nE(Ln  -  *  0  , 

n 

as  n  ■*  ®;  combining  (5.36)  and  (5.37),  we  get  (5.35i).  For  (ii),  we  use 
the  uniform  integrability  to  obtain 

2 

nE(s  -a)  +  k  ; 
n 

this  evidently  implies  that 

2  2  2  2 
n  E(L  -  EL  )  ♦  4z  (6)  k  , 

n  n 


proving  (ii). 

Me  conclude  that  the  end-point  variability  of  the  regenerative 
confidence  interval  is  of  order  n"1,  as  opposed  to  n~1/2  for 
standardized  time  series. 


i 


.mi 
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6 


SUMMARY 


Building  on  work  of  SCHRUBEN,  we  have  developed  a  general  framework 
for  Che  analysis  of  standardized  time  series.  Under  mild  assumptions  on 
the  output  process  (see  (3.1)),  the  method  of  standardized  time  series 
produces  asymptotically  valid  confidence  intevals  for  steady-state 
parameters.  However,  these  intervals  are  asymptotically  larger  (see 
(5.16))  and  more  variable  (see  (5.33))  then  those  steady-state  intervals 
obtained  by  a  method  which  consistently  estimates  the  appropriate  steady- 
state  variance  constant  (such  as  the  regenerateive  method).  In  this  sense 
standardized  time  series  confidence  Intervals  are  asymptotically  less 
desirable  then  those  constructed  by  a  consistent  estimation. 

These  results  do  not,  however,  preclude  the  possibility  that  standard¬ 
ized  times  series  may  be  superior  in  certain  small  sample  context;  this 
remains  an  area  for  future  work. 
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